#Results of Table 1 in the manuscript
#The R code files Univariate-CopulaModels1minus0.85.r, Univariate-CopulaModels2minus0.85.r, Univariate-CopulaModels3minus0.85.r and Univariate-CopulaModels4minus0.85.r 
#reproduce simulation results (% bias and RMSE) when fitting the 250 simulated datasets using the univariate model, Gaussian and 
#the flexible recursive bivariate model (detonated as Univariate, Gaussian C and Flexible C) with coefficient of the unobserved confounder set 
#to -0.85 in the treatment equation. 

#Univariate-CopulaModels1minus0.85.r allows one to reproduce simulation results for the logit, probit and cloglog link functions 
#(denoted as Logit-Logit, Probit-Probit and Cloglog-Cloglog) with normally distributed unmeasured confounder. 
#Similarly, Univariate-CopulaModels2minus0.85.r, Univariate-CopulaModels3minus0.85.r and Univariate-CopulaModels4minus0.85.r allow one to reproduce 
#simulation results for the logit, probit and cloglog link functions with Student�s t, Chi2 and uniform distributions 
#for the unobserved confounder, respectively.

#The four R code files, GM1minus0.85.r, GM2minus0.85.r, GM3minus0.85.r and GM4minus0.85.r, reproduce % bias and RMSE for the Matching when 
#setting the coefficient of the unobserved confounder to -0.85 in the treatment equation. 
#They allow one to reproduce simulation results for the logit, probit and cloglog link functions with normal, Student�s t, 
#Chi2 and uniform distributions for the unobserved confounder, respectively.

#When running Univariate-CopulaModels1minus0.85.r Univariate-CopulaModels2minus0.85.r, Univariate-CopulaModels3minus0.85.r, 
#Univariate-CopulaModels4minus0.85.r, GM1minus0.85.r, GM2minus0.85.r, GM3minus0.85.r and GM4minus0.85.r 

source("Univariate-CopulaModels1minus0.85.r")
source("Univariate-CopulaModels2minus0.85.r")
source("Univariate-CopulaModels3minus0.85.r")
source("Univariate-CopulaModels4minus0.85.r")
source("GM1minus0.85.r")
source("GM2minus0.85.r")
source("GM3minus0.85.r")
source("GM4minus0.85.r")


#replicates are stored in the following R workspaces: resminus0.85un1, resminus0.85un2, resminus0.85un3 and resminus0.85un4, resbetaminus0.85unob1, 
#resbetaminus0.85unob2, resbetaminus0.85unob3 and resbetaminus0.85unob4. 


#After having loaded all workspaces and stored the results using the following code


load("resminus0.85un1.RData")
results1 <- results
load("resminus0.85un2.RData")
results2 <- results
load("resminus0.85un3.RData")
results3 <- results
load("resminus0.85un4.RData")
results4 <- results
load("resbetaminus0.85unob1.RData")
results5 <- results
load("resbetaminus0.85unob2.RData")
results6 <- results
load("resbetaminus0.85unob3.RData")
results7 <- results
load("resbetaminus0.85unob4.RData")
results8 <- results

#results in Table 1 can be obtained using

bias1 <-round(results1$bias)
rmse1 <-round(sqrt(results1$mse),2)
bias2 <-round(results2$bias)
rmse2 <-round(sqrt(results2$mse),2)
bias3 <-round(results3$bias)
rmse3 <-round(sqrt(results3$mse),2)
bias4 <-round(results4$bias)
rmse4 <-round(sqrt(results4$mse),2)

bias5 <-round(results5$bias)
rmse5 <-round(sqrt(results5$mse),2)
bias6 <-round(results6$bias)
rmse6 <-round(sqrt(results6$mse),2)
bias7 <-round(results7$bias)
rmse7 <-round(sqrt(results7$mse),2)
bias8 <-round(results8$bias)
rmse8 <-round(sqrt(results8$mse),2)

bias <- matrix(NA, nrow=4, ncol=12)

bias[1,1] <- bias1[3]
bias[1,2] <- bias2[3]
bias[1,3] <- bias3[3]
bias[1,4] <- bias4[3]

bias[1,5] <- bias1[1]
bias[1,6] <- bias2[1]
bias[1,7] <- bias3[1]
bias[1,8] <- bias4[1]

bias[1,9] <- bias1[5]
bias[1,10] <- bias2[5]
bias[1,11] <- bias3[5]
bias[1,12] <- bias4[5]

bias[3,1] <- bias1[4]
bias[3,2] <- bias2[4]
bias[3,3] <- bias3[4]
bias[3,4] <- bias4[4]

bias[3,5] <- bias1[2]
bias[3,6] <- bias2[2]
bias[3,7] <- bias3[2]
bias[3,8] <- bias4[2]

bias[3,9] <- bias1[6]
bias[3,10] <- bias2[6]
bias[3,11] <- bias3[6]
bias[3,12] <- bias4[6]

bias[4,1] <- bias1[8]
bias[4,2] <- bias2[8]
bias[4,3] <- bias3[8]
bias[4,4] <- bias4[8]

bias[4,5] <- bias1[7]
bias[4,6] <- bias2[7]
bias[4,7] <- bias3[7]
bias[4,8] <- bias4[7]

bias[4,9] <- bias1[9]
bias[4,10] <- bias2[9]
bias[4,11] <- bias3[9]
bias[4,12] <- bias4[9]

bias[2,1] <- bias5[1]
bias[2,2] <- bias6[1]
bias[2,3] <- bias7[1]
bias[2,4] <- bias8[1]

bias[2,5] <- bias5[2]
bias[2,6] <- bias6[2]
bias[2,7] <- bias7[2]
bias[2,8] <- bias8[2]

bias[2,9] <- bias5[3]
bias[2,10] <- bias6[3]
bias[2,11] <- bias7[3]
bias[2,12] <- bias8[3]

rmse <- matrix(NA, nrow=4, ncol=12)

rmse[1,1] <- rmse1[3]
rmse[1,2] <- rmse2[3]
rmse[1,3] <- rmse3[3]
rmse[1,4] <- rmse4[3]

rmse[1,5] <- rmse1[1]
rmse[1,6] <- rmse2[1]
rmse[1,7] <- rmse3[1]
rmse[1,8] <- rmse4[1]

rmse[1,9] <- rmse1[5]
rmse[1,10] <- rmse2[5]
rmse[1,11] <- rmse3[5]
rmse[1,12] <- rmse4[5]

rmse[3,1] <- rmse1[4]
rmse[3,2] <- rmse2[4]
rmse[3,3] <- rmse3[4]
rmse[3,4] <- rmse4[4]

rmse[3,5] <- rmse1[2]
rmse[3,6] <- rmse2[2]
rmse[3,7] <- rmse3[2]
rmse[3,8] <- rmse4[2]

rmse[3,9] <- rmse1[6]
rmse[3,10] <- rmse2[6]
rmse[3,11] <- rmse3[6]
rmse[3,12] <- rmse4[6]

rmse[4,1] <- rmse1[8]
rmse[4,2] <- rmse2[8]
rmse[4,3] <- rmse3[8]
rmse[4,4] <- rmse4[8]

rmse[4,5] <- rmse1[7]
rmse[4,6] <- rmse2[7]
rmse[4,7] <- rmse3[7]
rmse[4,8] <- rmse4[7]

rmse[4,9] <- rmse1[9]
rmse[4,10] <- rmse2[9]
rmse[4,11] <- rmse3[9]
rmse[4,12] <- rmse4[9]

rmse[2,1] <- rmse5[1]
rmse[2,2] <- rmse6[1]
rmse[2,3] <- rmse7[1]
rmse[2,4] <- rmse8[1]

rmse[2,5] <- rmse5[2]
rmse[2,6] <- rmse6[2]
rmse[2,7] <- rmse7[2]
rmse[2,8] <- rmse8[2]

rmse[2,9] <- rmse5[3]
rmse[2,10] <- rmse6[3]
rmse[2,11] <- rmse7[3]
rmse[2,12] <- rmse8[3]


bias

rmse



#Results reported in Tables 1 and 2 of the Appendix can be obtained in a similar fashion.

source("Univariate-CopulaModels10.r")
source("Univariate-CopulaModels20.r")
source("Univariate-CopulaModels30.r")
source("Univariate-CopulaModels40.r")
source("GM10.r")
source("GM20.r")
source("GM30.r")
source("GM40.r")


source("Univariate-CopulaModels1minus1.5.r")
source("Univariate-CopulaModels2minus1.5.r")
source("Univariate-CopulaModels3minus1.5.r")
source("Univariate-CopulaModels4minus1.5.r")
source("GM1minus1.5.r")
source("GM2minus1.5.r")
source("GM3minus1.5.r")
source("GM4minus1.5.r")



#Table 1 of Appendix:

load("res0un1.RData")
results1 <- results
load("res0un2.RData")
results2 <- results
load("res0un3.RData")
results3 <- results
load("res0un4.RData")
results4 <- results
load("resbeta0unob1.RData")
results5 <- results
load("resbeta0unob2.RData")
results6 <- results
load("resbeta0unob3.RData")
results7 <- results
load("resbeta0unob4.RData")
results8 <- results
bias1 <-round(results1$bias)
rmse1 <-round(sqrt(results1$mse),2)
bias2 <-round(results2$bias)
rmse2 <-round(sqrt(results2$mse),2)
bias3 <-round(results3$bias)
rmse3 <-round(sqrt(results3$mse),2)
bias4 <-round(results4$bias)
rmse4 <-round(sqrt(results4$mse),2)

bias5 <-round(results5$bias)
rmse5 <-round(sqrt(results5$mse),2)
bias6 <-round(results6$bias)
rmse6 <-round(sqrt(results6$mse),2)
bias7 <-round(results7$bias)
rmse7 <-round(sqrt(results7$mse),2)
bias8 <-round(results8$bias)
rmse8 <-round(sqrt(results8$mse),2)

bias <- matrix(NA, nrow=4, ncol=12)

bias[1,1] <- bias1[3]
bias[1,2] <- bias2[3]
bias[1,3] <- bias3[3]
bias[1,4] <- bias4[3]

bias[1,5] <- bias1[1]
bias[1,6] <- bias2[1]
bias[1,7] <- bias3[1]
bias[1,8] <- bias4[1]

bias[1,9] <- bias1[5]
bias[1,10] <- bias2[5]
bias[1,11] <- bias3[5]
bias[1,12] <- bias4[5]

bias[3,1] <- bias1[4]
bias[3,2] <- bias2[4]
bias[3,3] <- bias3[4]
bias[3,4] <- bias4[4]

bias[3,5] <- bias1[2]
bias[3,6] <- bias2[2]
bias[3,7] <- bias3[2]
bias[3,8] <- bias4[2]

bias[3,9] <- bias1[6]
bias[3,10] <- bias2[6]
bias[3,11] <- bias3[6]
bias[3,12] <- bias4[6]

bias[4,1] <- bias1[8]
bias[4,2] <- bias2[8]
bias[4,3] <- bias3[8]
bias[4,4] <- bias4[8]

bias[4,5] <- bias1[7]
bias[4,6] <- bias2[7]
bias[4,7] <- bias3[7]
bias[4,8] <- bias4[7]

bias[4,9] <- bias1[9]
bias[4,10] <- bias2[9]
bias[4,11] <- bias3[9]
bias[4,12] <- bias4[9]

bias[2,1] <- bias5[1]
bias[2,2] <- bias6[1]
bias[2,3] <- bias7[1]
bias[2,4] <- bias8[1]

bias[2,5] <- bias5[2]
bias[2,6] <- bias6[2]
bias[2,7] <- bias7[2]
bias[2,8] <- bias8[2]

bias[2,9] <- bias5[3]
bias[2,10] <- bias6[3]
bias[2,11] <- bias7[3]
bias[2,12] <- bias8[3]

rmse <- matrix(NA, nrow=4, ncol=12)

rmse[1,1] <- rmse1[3]
rmse[1,2] <- rmse2[3]
rmse[1,3] <- rmse3[3]
rmse[1,4] <- rmse4[3]

rmse[1,5] <- rmse1[1]
rmse[1,6] <- rmse2[1]
rmse[1,7] <- rmse3[1]
rmse[1,8] <- rmse4[1]

rmse[1,9] <- rmse1[5]
rmse[1,10] <- rmse2[5]
rmse[1,11] <- rmse3[5]
rmse[1,12] <- rmse4[5]

rmse[3,1] <- rmse1[4]
rmse[3,2] <- rmse2[4]
rmse[3,3] <- rmse3[4]
rmse[3,4] <- rmse4[4]

rmse[3,5] <- rmse1[2]
rmse[3,6] <- rmse2[2]
rmse[3,7] <- rmse3[2]
rmse[3,8] <- rmse4[2]

rmse[3,9] <- rmse1[6]
rmse[3,10] <- rmse2[6]
rmse[3,11] <- rmse3[6]
rmse[3,12] <- rmse4[6]

rmse[4,1] <- rmse1[8]
rmse[4,2] <- rmse2[8]
rmse[4,3] <- rmse3[8]
rmse[4,4] <- rmse4[8]

rmse[4,5] <- rmse1[7]
rmse[4,6] <- rmse2[7]
rmse[4,7] <- rmse3[7]
rmse[4,8] <- rmse4[7]

rmse[4,9] <- rmse1[9]
rmse[4,10] <- rmse2[9]
rmse[4,11] <- rmse3[9]
rmse[4,12] <- rmse4[9]


rmse[2,1] <- rmse5[1]
rmse[2,2] <- rmse6[1]
rmse[2,3] <- rmse7[1]
rmse[2,4] <- rmse8[1]

rmse[2,5] <- rmse5[2]
rmse[2,6] <- rmse6[2]
rmse[2,7] <- rmse7[2]
rmse[2,8] <- rmse8[2]

rmse[2,9] <- rmse5[3]
rmse[2,10] <- rmse6[3]
rmse[2,11] <- rmse7[3]
rmse[2,12] <- rmse8[3]


bias 

rmse


#Table2 of Appendix:


load("resminus1.5un1.RData")
results1 <- results
load("resminus1.5un2.RData")
results2 <- results
load("resminus1.5un3.RData")
results3 <- results
load("resminus1.5un4.RData")
results4 <- results
load("resbetaminus1.5unob1.RData")
results5 <- results
load("resbetaminus1.5unob2.RData")
results6 <- results
load("resbetaminus1.5unob3.RData")
results7 <- results
load("resbetaminus1.5unob4.RData")
results8 <- results
bias1 <-round(results1$bias)
rmse1 <-round(sqrt(results1$mse),2)
bias2 <-round(results2$bias)
rmse2 <-round(sqrt(results2$mse),2)
bias3 <-round(results3$bias)
rmse3 <-round(sqrt(results3$mse),2)
bias4 <-round(results4$bias)
rmse4 <-round(sqrt(results4$mse),2)

bias5 <-round(results5$bias)
rmse5 <-round(sqrt(results5$mse),2)
bias6 <-round(results6$bias)
rmse6 <-round(sqrt(results6$mse),2)
bias7 <-round(results7$bias)
rmse7 <-round(sqrt(results7$mse),2)
bias8 <-round(results8$bias)
rmse8 <-round(sqrt(results8$mse),2)

bias <- matrix(NA, nrow=4, ncol=12)

bias[1,1] <- bias1[3]
bias[1,2] <- bias2[3]
bias[1,3] <- bias3[3]
bias[1,4] <- bias4[3]

bias[1,5] <- bias1[1]
bias[1,6] <- bias2[1]
bias[1,7] <- bias3[1]
bias[1,8] <- bias4[1]

bias[1,9] <- bias1[5]
bias[1,10] <- bias2[5]
bias[1,11] <- bias3[5]
bias[1,12] <- bias4[5]

bias[3,1] <- bias1[4]
bias[3,2] <- bias2[4]
bias[3,3] <- bias3[4]
bias[3,4] <- bias4[4]

bias[3,5] <- bias1[2]
bias[3,6] <- bias2[2]
bias[3,7] <- bias3[2]
bias[3,8] <- bias4[2]

bias[3,9] <- bias1[6]
bias[3,10] <- bias2[6]
bias[3,11] <- bias3[6]
bias[3,12] <- bias4[6]

bias[4,1] <- bias1[8]
bias[4,2] <- bias2[8]
bias[4,3] <- bias3[8]
bias[4,4] <- bias4[8]

bias[4,5] <- bias1[7]
bias[4,6] <- bias2[7]
bias[4,7] <- bias3[7]
bias[4,8] <- bias4[7]

bias[4,9] <- bias1[9]
bias[4,10] <- bias2[9]
bias[4,11] <- bias3[9]
bias[4,12] <- bias4[9]

bias[2,1] <- bias5[1]
bias[2,2] <- bias6[1]
bias[2,3] <- bias7[1]
bias[2,4] <- bias8[1]

bias[2,5] <- bias5[2]
bias[2,6] <- bias6[2]
bias[2,7] <- bias7[2]
bias[2,8] <- bias8[2]

bias[2,9] <- bias5[3]
bias[2,10] <- bias6[3]
bias[2,11] <- bias7[3]
bias[2,12] <- bias8[3]


rmse <- matrix(NA, nrow=4, ncol=12)

rmse[1,1] <- rmse1[3]
rmse[1,2] <- rmse2[3]
rmse[1,3] <- rmse3[3]
rmse[1,4] <- rmse4[3]

rmse[1,5] <- rmse1[1]
rmse[1,6] <- rmse2[1]
rmse[1,7] <- rmse3[1]
rmse[1,8] <- rmse4[1]

rmse[1,9] <- rmse1[5]
rmse[1,10] <- rmse2[5]
rmse[1,11] <- rmse3[5]
rmse[1,12] <- rmse4[5]

rmse[3,1] <- rmse1[4]
rmse[3,2] <- rmse2[4]
rmse[3,3] <- rmse3[4]
rmse[3,4] <- rmse4[4]

rmse[3,5] <- rmse1[2]
rmse[3,6] <- rmse2[2]
rmse[3,7] <- rmse3[2]
rmse[3,8] <- rmse4[2]

rmse[3,9] <- rmse1[6]
rmse[3,10] <- rmse2[6]
rmse[3,11] <- rmse3[6]
rmse[3,12] <- rmse4[6]

rmse[4,1] <- rmse1[8]
rmse[4,2] <- rmse2[8]
rmse[4,3] <- rmse3[8]
rmse[4,4] <- rmse4[8]

rmse[4,5] <- rmse1[7]
rmse[4,6] <- rmse2[7]
rmse[4,7] <- rmse3[7]
rmse[4,8] <- rmse4[7]

rmse[4,9] <- rmse1[9]
rmse[4,10] <- rmse2[9]
rmse[4,11] <- rmse3[9]
rmse[4,12] <- rmse4[9]

rmse[2,1] <- rmse5[1]
rmse[2,2] <- rmse6[1]
rmse[2,3] <- rmse7[1]
rmse[2,4] <- rmse8[1]

rmse[2,5] <- rmse5[2]
rmse[2,6] <- rmse6[2]
rmse[2,7] <- rmse7[2]
rmse[2,8] <- rmse8[2]

rmse[2,9] <- rmse5[3]
rmse[2,10] <- rmse6[3]
rmse[2,11] <- rmse7[3]
rmse[2,12] <- rmse8[3]

bias

rmse
